# Functions to support "Absolute and Relative Mobility: Two Frameworks for Connecting Intergenerational Mobility in Absolute and Relative Terms"
# (Sociological Methods and Research)

#####################
# NOTE TO USER:
####################

# This file is sourced in the file ``master_analysis.R''. Before running the master file, ensure that this function file and the data file are saved in user-set working directory

#####################
# General functions
####################

# Function to fix a bug in Hmisc::wtd.rank
wtd_rank <- function (x, weights = NULL, normwt = FALSE, na.rm = TRUE){
  if (!length(weights))
    return(rank(x, na.last = if (na.rm) NA else TRUE))
  tab <- Hmisc::wtd.table(x, weights, normwt = normwt, na.rm = na.rm)
  freqs <- tab$sum.of.weights
  r <- cumsum(freqs) - 0.5 * (freqs - 1)
  stats::approx(tab$x, r, xout = x, rule = 2)$y
}

# Function to calculate absolute mobility from relative mobility 
# (the parameter 'up' is the percentage increase in absolute income relative to parents, such that 20 = 20%)  
a.pop.fun <- function(mu.a,mu.p,sigma.a,sigma.p,rho,up) {
	mean.growth <- mu.a-mu.p  
	z.score.up0 <- mean.growth/sqrt(sigma.a^2+sigma.p^2-2*rho*sigma.a*sigma.p)
	up.to.log.diff <- log(up/100+1)
	z.score.up <- z.score.up0-up.to.log.diff
	a.mob <- 100*pnorm(z.score.up)
	out <- list(a.mob=a.mob) 
	return(out)
}

# Function to calculate absolute mobility from relative mobility,
# in different components  
a.pred.fun <- function(mu.a,mu.p,sigma.a,sigma.p,rho,x.p,n) {
	mean.growth <- mu.a-mu.p  
	mean.dist <- x.p-mu.p
	persist.ineq.dist <- rho*(sigma.a/sigma.p)*mean.dist 
	res.sd <- sqrt((1-rho^2)*sigma.a^2) 
	a.mob.expected <- mean.growth-mean.dist+persist.ineq.dist  
	a.mob.expected.sd.units <- a.mob.expected/res.sd
	mean.growth.sd.units <- mean.growth/res.sd
	mean.dist.sd.units <- mean.dist/res.sd
	persist.ineq.dist.sd.units <- persist.ineq.dist/res.sd 
	out <- list(mean.growth=mean.growth,mean.dist=mean.dist,
	            persist.ineq.dist=persist.ineq.dist,res.sd=res.sd, 
	            a.mob.expected=a.mob.expected,  
	            a.mob.expected.sd.units=a.mob.expected.sd.units,
	            mean.growth.sd.units=mean.growth.sd.units,mean.dist.sd.units=mean.dist.sd.units,
	            persist.ineq.dist.sd.units=persist.ineq.dist.sd.units) 
	return(out)
}

# Function to calculate absolute mobility from relative mobility,
# conditional on parental log income  
a.cond.fun <- function(mu.a,mu.p,sigma.a,sigma.p,rho,up,x.p) {
	mean.growth <- mu.a-mu.p 
	mean.dist.p <- x.p-mu.p 
	xp.above.mup <- x.p>mu.p
	xp.above.mua <- x.p>mu.a 
	a.mob.mat <- matrix(NA,nrow=length(x.p),ncol=length(rho)) 
	for(i in 1:(length(rho))) {
	  z.score.up0.num <- mean.growth-mean.dist.p+(mean.dist.p*rho[i]*sigma.a/sigma.p)
	  z.score.up0.denom <- sqrt((1-rho[i]^2)*sigma.a^2)
	  z.score.up0 <- z.score.up0.num/z.score.up0.denom 
	  up.to.log.diff <- log(up/100+1)
	  z.score.up <- z.score.up0-up.to.log.diff
	  a.mob <- 100*pnorm(z.score.up)
	  a.mob.mat[,i] <- a.mob 	    
	}
	out <- list(a.mob.mat=a.mob.mat,xp.above.mup=xp.above.mup,xp.above.mua=xp.above.mua) 
	return(out)
} 

# Function to calculate absolute mobility from relative mobility,
# allowing adult mean income to depend on relative mobility  
a.mua.pop.fun <- function(mu.a0,mu.p,sigma.a,sigma.p,rho,delta,up) {
	mean.growth.vec <- NULL
	mean.growth.pct.vec <- NULL
	a.mob.vec <- NULL
	for(i in 1:length(rho)) {
	mean.growth <- (mu.a0+delta*rho[i]*sigma.a/sigma.p)-mu.p 
	mean.growth.pct <- 100*(exp(mean.growth)-1) 
	z.score.up0 <- mean.growth/sqrt(sigma.a^2+sigma.p^2-2*rho[i]*sigma.a*sigma.p)
	up.to.log.diff <- log(up/100+1)
	z.score.up <- z.score.up0-up.to.log.diff
	a.mob <- 100*pnorm(z.score.up)
	a.mob.vec <- c(a.mob.vec,a.mob)
	mean.growth.vec <- c(mean.growth.vec,mean.growth)
	mean.growth.pct.vec <- c(mean.growth.pct.vec,mean.growth.pct)
	}
	out <- list(a.mob.vec=a.mob.vec,mean.growth.vec=mean.growth.vec,mean.growth.pct.vec=mean.growth.pct.vec) 
	return(out)
} 

# Function to simulate from multivariate t distribution  
t.sim.fun <- function(mu.a,mu.p,sigma.a,sigma.p,rho,n.add,df.add) {
  mu.vec <- c(mu.a,mu.p)
  pct.up <- NULL
  tsamp.col <- NULL
  for(i in 1:length(rho)) {
    sigma.mat <- matrix(c(sigma.a^2,rho[i]*sigma.a*sigma.p,rho[i]*sigma.a*sigma.p,sigma.p^2),ncol=2)
    tsamp <- rmvt(n.add,type="shifted",delta=mu.vec,sigma=sigma.mat*(df.add-2)/df.add,df=df.add)
    pct.up <- c(pct.up,100*sum(tsamp[,1]>tsamp[,2])/n.add)
    tsamp.col <- cbind(tsamp.col, tsamp)
  }
  out <- list(pct.up=pct.up,tsamp.col=tsamp.col)
  return(out)
}

# Function to generate random draws from multivariate gamma distribution using normal copula
# (see https://rdrr.io/rforge/lcmix/src/R/distributions.R and https://rdrr.io/rforge/lcmix/man/mvgamma.html)
rmvgamma <- function(n, shape=1, rate=1, corr=diag(length(shape))) {
	if(!is.matrix(corr) || !isSymmetric(corr))
		stop("'corr' must be a symmetric matrix")
	D = ncol(corr)
	Ds = length(shape)
	if(Ds > D)
		warning("'shape' longer than width of 'corr', truncating to fit")
	if(Ds != D)
		shape = rep(shape, length.out=D)
	Dr = length(rate)
	if(Dr > D)
		warning("'rate' longer than width of 'corr', truncating to fit")
	if(Dr != D)
		rate = rep(rate, length.out=D)
	if(D == 1) rgamma(n, shape, rate) 
	Z = rmvnorm(n, sigma=corr)
	cdf = pnorm(Z) 
	sapply(1:D, function(d) qgamma(cdf[,d], shape[d], rate[d]))
}

# Function to simulate from multivariate gamma distribution  
gamma.sim.fun <- function(mu.a,mu.p,sigma.a,sigma.p,rho,n.add,scale.add) {
  mu.vec <- c(mu.a,mu.p)
  sigma.vec <- c(sigma.a,sigma.p)
  shape.input <- ((mu.vec^2)/(sigma.vec^2))/scale.add
  rate.input <- mu.vec/(sigma.vec^2)
  pct.up <- NULL
  gsamp.col <- NULL
  for(i in 1:length(rho)) {
    rho.mat <- matrix(c(1,rho[i],rho[i],1),ncol=2)
    gsamp <- rmvgamma(n,shape=shape.input,rate=rate.input,corr=rho.mat)
    pct.up <- c(pct.up,100*sum(gsamp[,1]>gsamp[,2])/n.add)
    gsamp.col <- cbind(gsamp.col, gsamp)
  }
  out <- list(pct.up=pct.up,gsamp.col=gsamp.col)
}

####################
# Figure 1
####################

# Function to process some input for Figure 1
input.fig01.fun <- function(a.input) { 
	  out.tab <- cbind(c(wtd.mean(a.input$a.mob.expected[dat$origin_inc_above_mean==0],dat$weight[dat$origin_inc_above_mean==0]), 
                         wtd.mean(a.input$a.mob.expected.sd.units[dat$origin_inc_above_mean==0],dat$weight[dat$origin_inc_above_mean==0])),
                       c(wtd.mean(a.input$a.mob.expected[dat$origin_inc_above_mean==1],dat$weight[dat$origin_inc_above_mean==1]), 
                         wtd.mean(a.input$a.mob.expected.sd.units[dat$origin_inc_above_mean==1],dat$weight[dat$origin_inc_above_mean==1])))
	  colnames(out.tab) <- c("Below Parental Mean","Above Parental Mean")
	  rownames(out.tab) <- c("Expected, Log Units","Expected, SD units")
	  return(out.tab)
}

####################
# Figure 4
####################

# Function to create some input for Figure 4 (one group comparison)
fig04.fun.two.groups <- function(data.input,g1,g2) {
	  r.mob.obs <- data.input$dest_rank-data.input$origin_rank
	  a.mob.obs <- data.input$dest_inc_dol-data.input$origin_inc_dol
	  up.g1 <- wtd.mean(r.mob.obs[g1==1]>0 & a.mob.obs[g1==1]>0,data.input$weight[g1==1],na.rm=T)
	  down.g1 <- wtd.mean(r.mob.obs[g1==1]<0 & a.mob.obs[g1==1]<0,data.input$weight[g1==1],na.rm=T)
	  float.g1 <- wtd.mean(r.mob.obs[g1==1]<0 & a.mob.obs[g1==1]>0,data.input$weight[g1==1],na.rm=T)
	  up.g2 <- wtd.mean(r.mob.obs[g2==1]>0 & a.mob.obs[g2==1]>0,data.input$weight[g2==1],na.rm=T)
	  down.g2 <- wtd.mean(r.mob.obs[g2==1]<0 & a.mob.obs[g2==1]<0,data.input$weight[g2==1],na.rm=T)
	  float.g2 <- wtd.mean(r.mob.obs[g2==1]<0 & a.mob.obs[g2==1]>0,data.input$weight[g2==1],na.rm=T)                       
	  fig04.out <- cbind(up.g1-up.g2,down.g1-down.g2,float.g1-float.g2)
      colnames(fig04.out) <- c("Upward Diff.","Downward Diff.","Floating Diff.")  
      rownames(fig04.out) <- c("G1-G2")	
      return(fig04.out)
}

# Function to create some input for Figure 4 (all group comparisons)
fig04.fun.all <- function(data.input) {
	  fig04.wb <- fig04.fun.two.groups(data.input,data.input$white,data.input$black)
	  fig04.wh <- fig04.fun.two.groups(data.input,data.input$white,data.input$hispanic)
	  fig04.wo <- fig04.fun.two.groups(data.input,data.input$white,data.input$other)
	  fig04.bh <- fig04.fun.two.groups(data.input,data.input$black,data.input$hispanic)
	  fig04.bo <- fig04.fun.two.groups(data.input,data.input$black,data.input$other)
	  fig04.sn <- fig04.fun.two.groups(data.input,data.input$stable_par,1-data.input$stable_par)
	  fig04.nn <- fig04.fun.two.groups(data.input,data.input$native_par,1-data.input$native_par)
	  fig04.out.all <- 100*rbind(fig04.wb,fig04.wh,fig04.wo,fig04.bh,fig04.bo,fig04.sn,fig04.nn)
	  colnames(fig04.out.all) <- colnames(fig04.wb)
	  rownames(fig04.out.all) <- c("White - Black","White - Hispanic","White - Other","Black - Hispanic","Black - Other",
	                               "Stable Two Parents - Not","Native-Born Parents - Not")
	  return(fig04.out.all)                             
}

# Function to bootstrap some input for Figure 4 
fig04.fun.boot.pl <- function(data.input,above,nboot) {
  n.input <- dim(data.input)[1] 
  fig04.boot.list <- pbmclapply(1:nboot,function(i){ 
  	set.seed(i)
    samp.boot <- sample(1:n.input,n.input,replace=T)  
    data.input.boot <- data.input[samp.boot,]  
    mu.p.boot <- wtd.mean(data.input.boot$origin_inc_log,data.input.boot$weight)  
    data.input.boot <- data.input.boot %>%
         mutate(origin_rank = wtd_rank(origin_inc_log, weights = weight/1000, normwt = TRUE)/length(origin_inc_log),
                dest_rank = wtd_rank(dest_inc_log, weights = weight/1000, normwt = TRUE)/length(dest_inc_log),
                origin_inc_above_mean = origin_inc_log-mu.p.boot>0) 
	fig04.out.boot <- fig04.fun.all(data.input.boot[data.input.boot$origin_inc_above_mean==above,])               
    out <- list(fig04.out.boot=fig04.out.boot) 
    },mc.preschedule = FALSE, mc.set.seed = FALSE, mc.silent = FALSE, mc.cores = ncores)
  return(fig04.boot.list)
}

# Function to collect some bootstrapped output for Figure 4
fig04.fun.boot.collect <- function(fig04.boot.out,nboot) { 
	for(j in 1:3) {
	  assign(paste("col",j,sep="."),seq(j,3*nboot,3))
	}
	for(i in 1:1) {  
	  assign(paste("j.tab",i,sep="."),do.call(cbind,sapply(fig04.boot.out,"[",i)))
	  assign(paste("j.tab",i,"q025",sep="."),cbind(apply(get(paste("j.tab",i,sep="."))[,col.1],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.2],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.3],1,quantile,probs=.025,na.rm=T)))
	  assign(paste("j.tab",i,"q975",sep="."),cbind(apply(get(paste("j.tab",i,sep="."))[,col.1],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.2],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.3],1,quantile,probs=.975,na.rm=T))) 
	  tmp <- matrix(NA,nrow=7,ncol=3)
	  for(m in 1:7) {
	  for(n in 1:3) {
		tmp[m,n] <- paste("(",round(get(paste("j.tab",i,"q025",sep="."))[m,n],2),", ",
		                      round(get(paste("j.tab",i,"q975",sep="."))[m,n],2),")",sep="")
		}
		} 
	  assign(paste("j.tab",i,"ci",sep="."),tmp)                 
	} 
    out <- list(j.tab.1.q025=j.tab.1.q025,j.tab.1.q975=j.tab.1.q975,j.tab.1.ci=j.tab.1.ci) 
	return(out)            
}    

####################
# Table 1
####################

# Function to create Table 1 (point estimates)
tab01.fun <- function(a.input,origin.above.mean,weight) {
    tab01.out <- cbind(c(a.input$mean.growth,
                        -1*wtd.mean(a.input$mean.dist[origin.above.mean==0],weight[origin.above.mean==0]),
                        wtd.mean(a.input$persist.ineq.dist[origin.above.mean==0],weight[origin.above.mean==0]),
                        wtd.mean(a.input$a.mob.expected[origin.above.mean==0],weight[origin.above.mean==0])),
                       c(a.input$mean.growth, 
                         -1*wtd.mean(a.input$mean.dist[origin.above.mean==1],weight[origin.above.mean==1]),
                         wtd.mean(a.input$persist.ineq.dist[origin.above.mean==1],weight[origin.above.mean==1]),
                         wtd.mean(a.input$a.mob.expected[origin.above.mean==1],weight[origin.above.mean==1])),
                       c(a.input$mean.growth.sd.units,
                         -1*wtd.mean(a.input$mean.dist.sd.units[origin.above.mean==0],weight[origin.above.mean==0]),
                         wtd.mean(a.input$persist.ineq.dist.sd.units[origin.above.mean==0],weight[origin.above.mean==0]),
                         wtd.mean(a.input$a.mob.expected.sd.units[origin.above.mean==0],weight[origin.above.mean==0])),
                       c(a.input$mean.growth.sd.units, 
                         -1*wtd.mean(a.input$mean.dist.sd.units[origin.above.mean==1],weight[origin.above.mean==1]),
                         wtd.mean(a.input$persist.ineq.dist.sd.units[origin.above.mean==1],weight[origin.above.mean==1]),
                         wtd.mean(a.input$a.mob.expected.sd.units[origin.above.mean==1],weight[origin.above.mean==1])))
    colnames(tab01.out) <- c("Below Parental Mean","Above Parental Mean","Below Parental Mean (SD Units)","Above Parental Mean (SD Units)")
    rownames(tab01.out) <- c("Shared Growth","-Mean Reversion","Halted Reversion","Total")	
    return(tab01.out)	
}

# Function to bootstrap Table 1
tab01.fun.boot.pl <- function(data.input,nboot) {
  n.input <- dim(data.input)[1] 
  tab01.boot.list <- pbmclapply(1:nboot,function(i){ 
  	set.seed(i)
    samp.boot <- sample(1:n.input,n.input,replace=T)  
    data.input.boot <- data.input[samp.boot,]  
    mu.a.boot <- wtd.mean(data.input.boot$dest_inc_log,data.input.boot$weight)
	mu.p.boot <- wtd.mean(data.input.boot$origin_inc_log,data.input.boot$weight)
	sigma.a.boot <- sqrt(wtd.var(data.input.boot$dest_inc_log,data.input.boot$weight))
	sigma.p.boot <- sqrt(wtd.var(data.input.boot$origin_inc_log,data.input.boot$weight))
	rho.boot.reg <- coef(lm(I(scale(dest_inc_log))~I(scale(origin_inc_log)),weights=weight,data=data.input.boot))[2]  
    origin.above.mean.boot <- data.input.boot$origin_inc_log-mu.p.boot>0 
	a.pred.boot <- a.pred.fun(mu.a.boot,mu.p.boot,sigma.a.boot,sigma.p.boot, 
	                          rho.boot.reg,data.input.boot$origin_inc_log,n.input) 
    out.tab <- tab01.fun(a.pred.boot,origin.above.mean.boot,data.input.boot$weight)               
    out <- list(out.tab=out.tab)
    },mc.preschedule = FALSE, mc.set.seed = FALSE, mc.silent = FALSE, mc.cores = ncores)
  return(tab01.boot.list)
}

# Function to collect Table 1 bootstrapped output 
tab01.fun.boot.collect <- function(tab01.boot.output,nboot) { 
	for(j in 1:4) {
	  assign(paste("col",j,sep="."),seq(j,4*nboot,4))
	}
	for(i in 1:1) {  
	  assign(paste("j.tab",i,sep="."),do.call(cbind,sapply(tab01.boot.output,"[",i)))  
	  assign(paste("j.tab",i,"q025",sep="."),cbind(apply(get(paste("j.tab",i,sep="."))[,col.1],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.2],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.3],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.4],1,quantile,probs=.025,na.rm=T)))
	  assign(paste("j.tab",i,"q975",sep="."),cbind(apply(get(paste("j.tab",i,sep="."))[,col.1],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.2],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.3],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.4],1,quantile,probs=.975,na.rm=T))) 
	  tmp <- matrix(NA,nrow=4,ncol=4)
	  for(m in 1:4) {
	  for(n in 1:4) {
		tmp[m,n] <- paste("(",round(get(paste("j.tab",i,"q025",sep="."))[m,n],2),", ",
		                      round(get(paste("j.tab",i,"q975",sep="."))[m,n],2),")",sep="")
		}
		} 
	  assign(paste("j.tab",i,"ci",sep="."),tmp)                
	} 
    out <- list(j.tab.1.q025=j.tab.1.q025,j.tab.1.q975=j.tab.1.q975,j.tab.1.ci=j.tab.1.ci)  
	return(out)            
}  

####################
# Table 2
####################

# Function to create Table 2 (point estimates)
tab02.fun <- function(data.input) {
	  r.mob.obs <- data.input$dest_rank-data.input$origin_rank
	  a.mob.obs <- data.input$dest_inc_dol-data.input$origin_inc_dol
	  joint.quad <- ifelse(r.mob.obs>0 & a.mob.obs>0,1,
                           ifelse(r.mob.obs>0 & a.mob.obs<0,2,
                           ifelse(r.mob.obs<0 & a.mob.obs<0,3,4))) 
	  tab02.out <- rbind(c(100*wtd.mean(r.mob.obs>0 & a.mob.obs>0,data.input$weight,na.rm=T),
                              100*wtd.mean(r.mob.obs[joint.quad==1],data.input$weight[joint.quad==1]),
                              wtd.mean(a.mob.obs[joint.quad==1],data.input$weight[joint.quad==1]),        
                              100*wtd.mean(data.input$origin_rank[joint.quad==1],data.input$weight[joint.quad==1]), 
                              wtd.mean(data.input$origin_inc_dol[joint.quad==1],data.input$weight[joint.quad==1])),
                             c(100*wtd.mean(r.mob.obs<0 & a.mob.obs<0,data.input$weight,na.rm=T),
                               100*wtd.mean(r.mob.obs[joint.quad==3],data.input$weight[joint.quad==3]),
                               wtd.mean(a.mob.obs[joint.quad==3],data.input$weight[joint.quad==3]),        
                               100*wtd.mean(data.input$origin_rank[joint.quad==3],data.input$weight[joint.quad==3]), 
                               wtd.mean(data.input$origin_inc_dol[joint.quad==3],data.input$weight[joint.quad==3])), 
                             c(100*wtd.mean(r.mob.obs<0 & a.mob.obs>0,data.input$weight,na.rm=T),
                               100*wtd.mean(r.mob.obs[joint.quad==4],data.input$weight[joint.quad==4]),
                               wtd.mean(a.mob.obs[joint.quad==4],data.input$weight[joint.quad==4]),        
                               100*wtd.mean(data.input$origin_rank[joint.quad==4],data.input$weight[joint.quad==4]), 
                               wtd.mean(data.input$origin_inc_dol[joint.quad==4],data.input$weight[joint.quad==4])),
                             c(100*wtd.mean(r.mob.obs>0 & a.mob.obs<0,data.input$weight,na.rm=T),
                               100*wtd.mean(r.mob.obs[joint.quad==2],data.input$weight[joint.quad==2]),
                               wtd.mean(a.mob.obs[joint.quad==2],data.input$weight[joint.quad==2]),        
                               100*wtd.mean(data.input$origin_rank[joint.quad==2],data.input$weight[joint.quad==2]), 
                               wtd.mean(data.input$origin_inc_dol[joint.quad==2],data.input$weight[joint.quad==2])))
      colnames(tab02.out) <- c("%","Mean Rank $Delta$","Mean Abs. $Delta$","Mean Origin Rank","Mean Origin Inc.")  
      rownames(tab02.out) <- c("Upwardly Mobile","Downwardly Mobile","Floating with Rising Tide","Emerging from Retreating Tide")	
      return(tab02.out)
}

# Function to bootstrap Table 2 
tab02.fun.boot.pl <- function(data.input,nboot) {
  n.input <- dim(data.input)[1] 
  tab02.boot.list <- pbmclapply(1:nboot,function(i){ 
  	set.seed(i)
    samp.boot <- sample(1:n.input,n.input,replace=T)  
    data.input.boot <- data.input[samp.boot,] 
    data.input.boot <- data.input.boot %>%
         mutate(origin_rank = wtd_rank(origin_inc_log, weights = weight/1000, normwt = TRUE)/length(origin_inc_log),
                dest_rank = wtd_rank(dest_inc_log, weights = weight/1000, normwt = TRUE)/length(dest_inc_log)) 
	tab02.out.boot <- tab02.fun(data.input.boot)              
    out <- list(tab02.out.boot=tab02.out.boot)
    },mc.preschedule = FALSE, mc.set.seed = FALSE, mc.silent = FALSE, mc.cores = ncores)
  return(tab02.boot.list)
}

# Function to collect Table 2 bootstrapped output 
tab02.fun.boot.collect <- function(tab02.boot.out,nboot) { 
	for(j in 1:5) {
	  assign(paste("col",j,sep="."),seq(j,5*nboot,5))
	}
	for(i in 1:1) {  
	  assign(paste("j.tab",i,sep="."),do.call(cbind,sapply(tab02.boot.out,"[",i)))
	  assign(paste("j.tab",i,"q025",sep="."),cbind(apply(get(paste("j.tab",i,sep="."))[,col.1],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.2],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.3],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.4],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.5],1,quantile,probs=.025,na.rm=T)))
	  assign(paste("j.tab",i,"q975",sep="."),cbind(apply(get(paste("j.tab",i,sep="."))[,col.1],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.2],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.3],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.4],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.5],1,quantile,probs=.975,na.rm=T))) 
	  tmp <- matrix(NA,nrow=4,ncol=5)
	  for(m in 1:4) {
	  for(n in 1:5) {
		tmp[m,n] <- paste("(",round(get(paste("j.tab",i,"q025",sep="."))[m,n],2),", ",
		                      round(get(paste("j.tab",i,"q975",sep="."))[m,n],2),")",sep="")
		}
		} 
	  assign(paste("j.tab",i,"ci",sep="."),tmp)                 
	} 
    out <- list(j.tab.1.q025=j.tab.1.q025,j.tab.1.q975=j.tab.1.q975,j.tab.1.ci=j.tab.1.ci) 
	return(out)            
}    

# Function to combine Table 2 bootstrapped output with point estimates
tab02.fun.boot.combine <- function(tab02.point,tab02.boot.coll.out,ntable) { 
	tmp <- tab02.boot.coll.out[[3*ntable]]
	tab02.point.ci <- t(data.frame(round(tab02.point[1,],2),tmp[1,],
	                         round(tab02.point[2,],2),tmp[2,],
	                         round(tab02.point[3,],2),tmp[3,],
	                         round(tab02.point[4,],2),tmp[4,]))
	rownames(tab02.point.ci) <- c("Upwardly Mobile","CI1","Downwardly Mobile","CI2",
	                              "Floating with Rising Tide","CI3","Emerging from Retreating Tide","CI4") 
	return(tab02.point.ci)
}

####################
# Table 3
####################

# Function to create Table 3 (comparison between two groups, point estimates)
tab03.fun.two.groups <- function(data.input,g1,g2) {
	  r.mob.obs <- data.input$dest_rank-data.input$origin_rank
	  a.mob.obs <- data.input$dest_inc_dol-data.input$origin_inc_dol
	  down.g1 <- wtd.mean(r.mob.obs[g1==1]<0 & a.mob.obs[g1==1]<0,data.input$weight[g1==1],na.rm=T)
	  float.g1 <- wtd.mean(r.mob.obs[g1==1]<0 & a.mob.obs[g1==1]>0,data.input$weight[g1==1],na.rm=T)
	  float.share.g1 <- 100*float.g1/(float.g1+down.g1)
	  down.g2 <- wtd.mean(r.mob.obs[g2==1]<0 & a.mob.obs[g2==1]<0,data.input$weight[g2==1],na.rm=T)
	  float.g2 <- wtd.mean(r.mob.obs[g2==1]<0 & a.mob.obs[g2==1]>0,data.input$weight[g2==1],na.rm=T)
	  float.share.g2 <- 100*float.g2/(float.g2+down.g2)                        
	  tab03.out <- cbind(float.share.g1,float.share.g2,float.share.g1-float.share.g2,100*float.share.g1/float.share.g2-100)
      colnames(tab03.out) <- c("G1 %","G2 %","Diff.","Ratio")  
      rownames(tab03.out) <- c("G1-G2")	
      return(tab03.out)
}

# Function to create Table 3 (all group comparisons, point estimates)
tab03.fun.all <- function(data.input) {
	  tab03.wb <- tab03.fun.two.groups(data.input,data.input$white,data.input$black)
	  tab03.wh <- tab03.fun.two.groups(data.input,data.input$white,data.input$hispanic)
	  tab03.wo <- tab03.fun.two.groups(data.input,data.input$white,data.input$other)
	  tab03.bh <- tab03.fun.two.groups(data.input,data.input$black,data.input$hispanic)
	  tab03.bo <- tab03.fun.two.groups(data.input,data.input$black,data.input$other)
	  tab03.sn <- tab03.fun.two.groups(data.input,data.input$stable_par,1-data.input$stable_par)
	  tab03.nn <- tab03.fun.two.groups(data.input,data.input$native_par,1-data.input$native_par)
	  tab03.out.all <- rbind(tab03.wb,tab03.wh,tab03.wo,tab03.bh,tab03.bo,tab03.sn,tab03.nn)
	  colnames(tab03.out.all) <- colnames(tab03.wb)
	  rownames(tab03.out.all) <- c("White - Black","White - Hispanic","White - Other","Black - Hispanic","Black - Other",
	                               "Stable Two Parents - Not","Native-Born Parents - Not")
	  return(tab03.out.all)                             
}

# Function to bootstrap Table 3 
tab03.fun.boot.pl <- function(data.input,nboot) {
  n.input <- dim(data.input)[1] 
  tab03.boot.list <- pbmclapply(1:nboot,function(i){ 
  	set.seed(i)
    samp.boot <- sample(1:n.input,n.input,replace=T)  
    data.input.boot <- data.input[samp.boot,] 
    data.input.boot <- data.input.boot %>%
         mutate(origin_rank = wtd_rank(origin_inc_log, weights = weight/1000, normwt = TRUE)/length(origin_inc_log),
                dest_rank = wtd_rank(dest_inc_log, weights = weight/1000, normwt = TRUE)/length(dest_inc_log)) 
	tab03.out.boot <- tab03.fun.all(data.input.boot)              
    out <- list(tab03.out.boot=tab03.out.boot)
    },mc.preschedule = FALSE, mc.set.seed = FALSE, mc.silent = FALSE, mc.cores = ncores)
  return(tab03.boot.list)
}

# Function to collect Table 3 bootstrapped output 
tab03.fun.boot.collect <- function(tab03.boot.out,nboot) { 
	for(j in 1:4) {
	  assign(paste("col",j,sep="."),seq(j,4*nboot,4))
	}
	for(i in 1:1) {  
	  assign(paste("j.tab",i,sep="."),do.call(cbind,sapply(tab03.boot.out,"[",i)))
	  assign(paste("j.tab",i,"q025",sep="."),cbind(apply(get(paste("j.tab",i,sep="."))[,col.1],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.2],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.3],1,quantile,probs=.025,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.4],1,quantile,probs=.025,na.rm=T)))
	  assign(paste("j.tab",i,"q975",sep="."),cbind(apply(get(paste("j.tab",i,sep="."))[,col.1],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.2],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.3],1,quantile,probs=.975,na.rm=T),
	                                            apply(get(paste("j.tab",i,sep="."))[,col.4],1,quantile,probs=.975,na.rm=T))) 
	  tmp <- matrix(NA,nrow=7,ncol=4)
	  for(m in 1:7) {
	  for(n in 1:4) {
		tmp[m,n] <- paste("(",round(get(paste("j.tab",i,"q025",sep="."))[m,n],2),", ",
		                      round(get(paste("j.tab",i,"q975",sep="."))[m,n],2),")",sep="")
		}
		} 
	  assign(paste("j.tab",i,"ci",sep="."),tmp)                 
	} 
    out <- list(j.tab.1.q025=j.tab.1.q025,j.tab.1.q975=j.tab.1.q975,j.tab.1.ci=j.tab.1.ci) 
	return(out)            
}    
